## Do perceptions of urban bias in solar farms drive rural opposition to solar energy? 
## A survey experiment in South Korea

## Author: Inhwan Ko (University of Nevada, Reno)
## Last Updated on April 29, 202

pacman::p_load(MASS, readxl, tidyverse, car, simcf, tile, lmtest, sandwich,
               rstan, rstanarm, bayesplot, HDInterval)

data <- read_excel("data.xlsx")

data$gender <- recode(data$gender, "1=1;2=0")
data$solar <- as.numeric(data$solar)
data$nucl <- as.numeric(data$nucl)
data$q1_nuclecon <- as.numeric(data$q1_nuclecon)
data$q2_polbias <- as.numeric(data$q2_polbias)
data$q3_gap <- as.numeric(data$q3_gap)
data$gender <- as.numeric(data$gender)
data$income <- as.numeric(data$income)
data$educ <- as.numeric(data$educ)
data$ideo <- as.numeric(data$ideo)
data$age <- as.numeric(data$age)
data$age <- 2022-data$age

data$treat <- NA
data$treat[data$treatA==1] <- "A"
data$treat[data$treatB==1] <- "B"
data$treat[data$treatC==1] <- "C"

efftdata <- data %>% filter(efft==1) # filtering effective responses only

town1 <- efftdata %>% filter(area==1)
town2 <- efftdata %>% filter(area==2)

###############################################################################
################################# OLS Results #################################
###############################################################################

lm1 <- lm(solar ~ treatB + treatC +
               q1_nuclecon + q2_polbias + q3_gap + ideo + educ + age + gender + income, 
             efftdata)
summary(lm1)
round(coeftest(lm1, vcov=vcovHC(lm1, type="HC4")),3)

lm2 <- lm(solar ~ treatB + treatC +
            q1_nuclecon + q2_polbias + q3_gap + ideo + educ + age + gender + income, 
             town1)
summary(lm2)
round(coeftest(lm2, vcov=vcovHC(lm2, type="HC4")), 3)

lm3 <- lm(solar ~ treatB + treatC +
            q1_nuclecon + q2_polbias + q3_gap + ideo + educ + age + gender + income, 
          town2)
summary(lm3)
round(coeftest(lm3, vcov=vcovHC(lm3, type="HC4")), 3)

## OLS without potential post-treatment bias-inducing variables

lm1_pt <- lm(solar ~ treatB + treatC + q1_nuclecon + gender + income + educ + ideo + age, 
          efftdata)
summary(lm1_pt)
round(coeftest(lm1_pt, vcov=vcovHC(lm1_pt, type="HC4")),3)

lm2_pt <- lm(solar ~ treatB + treatC + q1_nuclecon + gender + income + educ + ideo + age, 
     town1)
summary(lm2_pt)
round(coeftest(lm2_pt, vcov=vcovHC(lm2_pt, type="HC4")), 3)

lm3_pt <- lm(solar ~ treatB + treatC + q1_nuclecon + gender + income + educ + ideo + age, 
     town2)
summary(lm3_pt)
round(coeftest(lm3_pt, vcov=vcovHC(lm3_pt, type="HC4")), 3)

###############################################################################
################ OLS with interaction terms with Bayesian approach ############
###############################################################################

rstan_options(iter=5000)
rstan_options(chains=4)

### Original results ###
bayesdata1 <- efftdata
colnames(bayesdata1)[4:5] <- c("Vignette 1", "Vignette 2")
colnames(bayesdata1)[12:19] <- c("Nuclear Reliance Question", "Urban Bias Question",
                                 "Rural-Urban Gap Question", "Gender", "Household Income",
                                 "Education", "Political Ideology", "Age")

formula1 <- solar ~ `Vignette 1` + `Vignette 2` +
  `Nuclear Reliance Question` + `Urban Bias Question` + `Rural-Urban Gap Question` + 
  `Political Ideology` + `Education` + `Age` + `Gender` + `Household Income`

prior_mean <- lm1$coefficients[2:11]
prior_scale <- sqrt(diag(vcovHC(lm1, type="HC4")))[2:11]
intercept_mean <- lm1$coefficients[1]
intercept_scale <- sqrt(diag(vcovHC(lm1, type="HC4")))[1]

stan_lm1 <- stan_glm(formula1, family=gaussian(), data=bayesdata1,
                     prior=normal(
                       prior_mean, prior_scale, autoscale=F),
                     prior_intercept=normal(
                       intercept_mean, intercept_scale, autoscale=F))

summary(stan_lm1, probs=c(0.025,0.975))
plot(stan_lm1)
mcmc_trace(stan_lm1, pars=c("Age","sigma"))

bayesdata2 <- bayesdata1 %>% filter(area==1)
bayesdata3 <- bayesdata1 %>% filter(area==2)

prior_mean <- lm2$coefficients[2:11]
prior_scale <- sqrt(diag(vcovHC(lm2, type="HC4")))[2:11]
intercept_mean <- lm2$coefficients[1]
intercept_scale <- sqrt(diag(vcovHC(lm2, type="HC4")))[1]

stan_lm2 <- stan_glm(formula1, family=gaussian(), data=bayesdata2,
                     prior=normal(
                       prior_mean, prior_scale, autoscale=F),
                     prior_intercept=normal(
                       intercept_mean, intercept_scale, autoscale=F))
summary(stan_lm2, probs=c(0.025,0.975))
plot(stan_lm2)
mcmc_trace(stan_lm2)

prior_mean <- lm3$coefficients[2:11]
prior_scale <- sqrt(diag(vcovHC(lm3, type="HC4")))[2:11]
intercept_mean <- lm3$coefficients[1]
intercept_scale <- sqrt(diag(vcovHC(lm3, type="HC4")))[1]

stan_lm3 <- stan_glm(formula1, family=gaussian(), data=bayesdata3,
                     prior=normal(
                       prior_mean, prior_scale, autoscale=F),
                     prior_intercept=normal(
                       intercept_mean, intercept_scale, autoscale=F))
summary(stan_lm3, probs=c(0.025,0.975))
plot(stan_lm3)
mcmc_trace(stan_lm3)

## interaction terms
formula2 <- solar ~ `Vignette 1` + `Vignette 2` +
  `Nuclear Reliance Question` + `Urban Bias Question` + `Rural-Urban Gap Question` + 
  `Political Ideology` + `Education` + `Age` + `Gender` + `Household Income` +
  `Vignette 1`*`Urban Bias Question` + `Vignette 2`*`Urban Bias Question`

stan_lm1b <- stan_glm(formula2, family=gaussian(), data=bayesdata1)
stan_lm2b <- stan_glm(formula2, family=gaussian(), data=bayesdata2)
stan_lm3b <- stan_glm(formula2, family=gaussian(), data=bayesdata3)

summary(stan_lm1b, probs=c(0.025,0.975))
summary(stan_lm2b, probs=c(0.025,0.975))
summary(stan_lm3b, probs=c(0.025,0.975))

formula3 <- solar ~ `Vignette 1` + `Vignette 2` +
  `Nuclear Reliance Question` + `Urban Bias Question` + `Rural-Urban Gap Question` + 
  `Political Ideology` + `Education` + `Age` + `Gender` + `Household Income` +
  `Vignette 1`*`Rural-Urban Gap Question` + `Vignette 2`*`Rural-Urban Gap Question`

stan_lm1c <- stan_glm(formula3, family=gaussian(), data=bayesdata1)
stan_lm2c <- stan_glm(formula3, family=gaussian(), data=bayesdata2)
stan_lm3c <- stan_glm(formula3, family=gaussian(), data=bayesdata3)

summary(stan_lm1c, probs=c(0.025,0.975))
summary(stan_lm2c, probs=c(0.025,0.975))
summary(stan_lm3c, probs=c(0.025,0.975))

###############################################################################
############################ Descriptive Statistics ###########################
###############################################################################

formula2 <- solar ~ `Vignette 1`*`Urban Bias Question` + 
  `Vignette 2`*`Urban Bias Question` +
  `Nuclear Reliance Question` + `Urban Bias Question` + `Rural-Urban Gap Question` + 
  `Political Ideology` + `Education` + `Age` + `Gender` + `Household Income`

stan_lm1a <- stan_glm(formula2, family=gaussian(), data=bayesdata1)
summary(stan_lm1a)
plot(stan_lm1a)
mcmc_trace(stan_lm1a)

stan_lm2a <- stan_glm(formula2, family=gaussian(), data=bayesdata2)
summary(stan_lm2a)
plot(stan_lm2a)
mcmc_trace(stan_lm2a)

stan_lm3a <- stan_glm(formula2, family=gaussian(), data=bayesdata3)
summary(stan_lm3a)
plot(stan_lm3a)
mcmc_trace(stan_lm3a)

### Interaction 2 ###
formula3 <- solar ~ `Vignette 1`*`Rural-Urban Gap Question` + 
  `Vignette 2`*`Rural-Urban Gap Question` +
  `Nuclear Reliance Question` + `Urban Bias Question` + `Rural-Urban Gap Question` + 
  `Political Ideology` + `Education` + `Age` + `Gender` + `Household Income`

stan_lm1b <- stan_glm(formula3, family=gaussian(), data=bayesdata1)
summary(stan_lm1b)
plot(stan_lm1b)
mcmc_trace(stan_lm1b)

stan_lm2b <- stan_glm(formula3, family=gaussian(), data=bayesdata2)
summary(stan_lm2b)
plot(stan_lm2b)
mcmc_trace(stan_lm2b)

stan_lm2c <- stan_glm(formula3, family=gaussian(), data=bayesdata3)
summary(stan_lm2c)
plot(stan_lm2c)
mcmc_trace(stan_lm2c)

#### Descriptive statistics ####
# response counts
data %>% filter(treatA==1 & area==1) %>% nrow()
data %>% filter(treatB==1 & area==1) %>% nrow()
data %>% filter(treatC==1 & area==1) %>% nrow()

data %>% filter(treatA==1 & area==2) %>% nrow()
data %>% filter(treatB==1 & area==2) %>% nrow()
data %>% filter(treatC==1 & area==2) %>% nrow()

data %>% filter(treatA==1 & area==1 & efft==0) %>% nrow()
data %>% filter(treatB==1 & area==1 & efft==0) %>% nrow()
data %>% filter(treatC==1 & area==1 & efft==0) %>% nrow()

data %>% filter(treatA==1 & area==2 & efft==0) %>% nrow()
data %>% filter(treatB==1 & area==2 & efft==0) %>% nrow()
data %>% filter(treatC==1 & area==2 & efft==0) %>% nrow()

data %>% filter(treatA==1 & area==1 & attcorrect==0) %>% nrow()
data %>% filter(treatB==1 & area==1 & attcorrect==0) %>% nrow()
data %>% filter(treatC==1 & area==1 & attcorrect==0) %>% nrow()

data %>% filter(treatA==1 & area==2 & attcorrect==0) %>% nrow()
data %>% filter(treatB==1 & area==2 & attcorrect==0) %>% nrow()
data %>% filter(treatC==1 & area==2 & attcorrect==0) %>% nrow()

# balance statistics
## town 1
town1 %>% filter(treatA==1) %>% select(q1_nuclecon) %>% 
  summarize(mean=mean(q1_nuclecon), sd=sd(q1_nuclecon))
town1 %>% filter(treatB==1) %>% select(q1_nuclecon) %>% 
  summarize(mean=mean(q1_nuclecon), sd=sd(q1_nuclecon))
town1 %>% filter(treatC==1) %>% select(q1_nuclecon) %>%
  summarize(mean=mean(q1_nuclecon), sd=sd(q1_nuclecon))
aov(q1_nuclecon ~ treat, town1) %>% summary()

town1 %>% filter(treatA==1) %>% select(q2_polbias) %>% 
  summarize(mean=mean(q2_polbias), sd=sd(q2_polbias))
town1 %>% filter(treatB==1) %>% select(q2_polbias) %>% 
  summarize(mean=mean(q2_polbias), sd=sd(q2_polbias))
town1 %>% filter(treatC==1) %>% select(q2_polbias) %>% 
  summarize(mean=mean(q2_polbias), sd=sd(q2_polbias))
aov(q2_polbias ~ treat, town1) %>% summary()

town1 %>% filter(treatA==1) %>% select(q3_gap) %>% 
  summarize(mean=mean(q3_gap), sd=sd(q3_gap))
town1 %>% filter(treatB==1) %>% select(q3_gap) %>%
  summarize(mean=mean(q3_gap), sd=sd(q3_gap))
town1 %>% filter(treatC==1) %>% select(q3_gap) %>% 
  summarize(mean=mean(q3_gap), sd=sd(q3_gap))
aov(q3_gap ~ treat, town1) %>% summary()

town1 %>% filter(treatA==1) %>% select(age) %>% 
  summarize(mean=mean(age), sd=sd(age))
town1 %>% filter(treatB==1) %>% select(age) %>% 
  summarize(mean=mean(age), sd=sd(age))
town1 %>% filter(treatC==1) %>% select(age) %>% 
  summarize(mean=mean(age), sd=sd(age))
aov(age ~ treat, town1) %>% summary()

town1 %>% filter(treatA==1) %>% select(gender) %>% 
  summarize(mean=mean(gender), sd=sd(gender))
town1 %>% filter(treatB==1) %>% select(gender) %>% 
  summarize(mean=mean(gender), sd=sd(gender))
town1 %>% filter(treatC==1) %>% select(gender) %>% 
  summarize(mean=mean(gender), sd=sd(gender))
aov(gender ~ treat, town1) %>% summary()

town1 %>% filter(treatA==1) %>% select(income) %>% 
  summarize(mean=mean(income), sd=sd(income))
town1 %>% filter(treatB==1) %>% select(income) %>% 
  summarize(mean=mean(income), sd=sd(income))
town1 %>% filter(treatC==1) %>% select(income) %>% 
  summarize(mean=mean(income), sd=sd(income))
aov(income ~ treat, town1) %>% summary()

town1 %>% filter(treatA==1) %>% select(educ) %>% 
  summarize(mean=mean(educ), sd=sd(educ))
town1 %>% filter(treatB==1) %>% select(educ) %>% 
  summarize(mean=mean(educ), sd=sd(educ))
town1 %>% filter(treatC==1) %>% select(educ) %>% 
  summarize(mean=mean(educ), sd=sd(educ))
aov(educ ~ treat, town1) %>% summary()

town1 %>% filter(treatA==1) %>% select(ideo) %>% 
  summarize(mean=mean(ideo), sd=sd(ideo))
town1 %>% filter(treatB==1) %>% select(ideo) %>% 
  summarize(mean=mean(ideo), sd=sd(ideo))
town1 %>% filter(treatC==1) %>% select(ideo) %>% 
  summarize(mean=mean(ideo), sd=sd(ideo))
aov(ideo ~ treat, town1) %>% summary()

## town 2
town2 %>% filter(treatA==1) %>% select(q1_nuclecon) %>% 
  summarize(mean=mean(q1_nuclecon), sd=sd(q1_nuclecon))
town2 %>% filter(treatB==1) %>% select(q1_nuclecon) %>% 
  summarize(mean=mean(q1_nuclecon), sd=sd(q1_nuclecon))
town2 %>% filter(treatC==1) %>% select(q1_nuclecon) %>%
  summarize(mean=mean(q1_nuclecon), sd=sd(q1_nuclecon))
aov(q1_nuclecon ~ treat, town2) %>% summary()

town2 %>% filter(treatA==1) %>% select(q2_polbias) %>% 
  summarize(mean=mean(q2_polbias), sd=sd(q2_polbias))
town2 %>% filter(treatB==1) %>% select(q2_polbias) %>% 
  summarize(mean=mean(q2_polbias), sd=sd(q2_polbias))
town2 %>% filter(treatC==1) %>% select(q2_polbias) %>% 
  summarize(mean=mean(q2_polbias), sd=sd(q2_polbias))
aov(q2_polbias ~ treat, town2) %>% summary()

town2 %>% filter(treatA==1) %>% select(q3_gap) %>% 
  summarize(mean=mean(q3_gap), sd=sd(q3_gap))
town2 %>% filter(treatB==1) %>% select(q3_gap) %>%
  summarize(mean=mean(q3_gap), sd=sd(q3_gap))
town2 %>% filter(treatC==1) %>% select(q3_gap) %>% 
  summarize(mean=mean(q3_gap), sd=sd(q3_gap))
aov(q3_gap ~ treat, town2) %>% summary()

town2 %>% filter(treatA==1) %>% select(age) %>% 
  summarize(mean=mean(age), sd=sd(age))
town2 %>% filter(treatB==1) %>% select(age) %>% 
  summarize(mean=mean(age), sd=sd(age))
town2 %>% filter(treatC==1) %>% select(age) %>% 
  summarize(mean=mean(age), sd=sd(age))
aov(age ~ treat, town2) %>% summary()

town2 %>% filter(treatA==1) %>% select(gender) %>% 
  summarize(mean=mean(gender), sd=sd(gender))
town2 %>% filter(treatB==1) %>% select(gender) %>% 
  summarize(mean=mean(gender), sd=sd(gender))
town2 %>% filter(treatC==1) %>% select(gender) %>% 
  summarize(mean=mean(gender), sd=sd(gender))
aov(gender ~ treat, town2) %>% summary()

town2 %>% filter(treatA==1) %>% select(income) %>% 
  summarize(mean=mean(income), sd=sd(income))
town2 %>% filter(treatB==1) %>% select(income) %>% 
  summarize(mean=mean(income), sd=sd(income))
town2 %>% filter(treatC==1) %>% select(income) %>% 
  summarize(mean=mean(income), sd=sd(income))
aov(income ~ treat, town2) %>% summary()

town2 %>% filter(treatA==1) %>% select(educ) %>% 
  summarize(mean=mean(educ), sd=sd(educ))
town2 %>% filter(treatB==1) %>% select(educ) %>% 
  summarize(mean=mean(educ), sd=sd(educ))
town2 %>% filter(treatC==1) %>% select(educ) %>% 
  summarize(mean=mean(educ), sd=sd(educ))
aov(educ ~ treat, town2) %>% summary()

town2 %>% filter(treatA==1) %>% select(ideo) %>% 
  summarize(mean=mean(ideo), sd=sd(ideo))
town2 %>% filter(treatB==1) %>% select(ideo) %>% 
  summarize(mean=mean(ideo), sd=sd(ideo))
town2 %>% filter(treatC==1) %>% select(ideo) %>% 
  summarize(mean=mean(ideo), sd=sd(ideo))
aov(ideo ~ treat, town2) %>% summary()


aov(q1_nuclecon ~ treat, efftdata) %>% summary()
aov(q2_polbias ~ treat, efftdata) %>% summary()
aov(q5_gap ~ treat, efftdata) %>% summary()
aov(age ~ treat, efftdata) %>% summary()
aov(gender ~ treat, efftdata) %>% summary()
aov(income ~ treat, efftdata) %>% summary()
aov(educ ~ treat, efftdata) %>% summary()
aov(ideo ~ treat, efftdata) %>% summary()

# balance check between attention-check failures and non-failures
data$attcorrect <- as.numeric(data$attcorrect)
aov(attcorrect ~ treat, data) %>% summary()

aov(q1_impact ~ attcorrect, data) %>% summary() # difference
mean(data$q1_impact[data$attcorrect==1], na.rm=T)
mean(data$q1_impact[data$attcorrect==0], na.rm=T)
t.test(data$q1_impact[data$attcorrect==1], data$q1_impact[data$attcorrect==0])

aov(q2_renewair ~ attcorrect, data) %>% summary()
aov(q3_nuclecon ~ attcorrect, data) %>% summary()

aov(q4_polbias ~ attcorrect, data) %>% summary() #
mean(data$q4_polbias[data$attcorrect==1], na.rm=T)
mean(data$q4_polbias[data$attcorrect==0], na.rm=T)
t.test(data$q4_polbias[data$attcorrect==1], data$q4_polbias[data$attcorrect==0])

aov(q5_gap ~ attcorrect, data) %>% summary()
aov(age ~ attcorrect, data) %>% summary()
aov(gender ~ attcorrect, data) %>% summary()
aov(income ~ attcorrect, data) %>% summary()
aov(educ ~ attcorrect, data) %>% summary()
aov(ideo ~ attcorrect, data) %>% summary() 
aov(reside ~ attcorrect, data) %>% summary()
aov(bornhere ~ attcorrect, data) %>% summary()
aov(relatseoul ~ attcorrect, data) %>% summary()

## map

library(sf)
firstmap <- st_read("ctprvn.shp")
secondmap <- st_read("sig.shp")
  
library(ggplot2)
secondmap.data <- as.data.frame(secondmap) 
Yeoju <- secondmap.data[115,]
Wando <- secondmap.data[183,]

firstmap[18,] <- Yeoju
firstmap[19,] <- Wando

firstmap$color <- ifelse(firstmap$CTPRVN_CD==11 |
                           firstmap$CTPRVN_CD==41670 |
                           firstmap$CTPRVN_CD==46890, 1, 0)

library(tmap)
map <- firstmap %>% 
  tm_shape() +
  tm_borders(lwd=0.5, alpha=1) +
  tm_fill("color", palette=c("white", "#5F190F"), style="cat",
          midpoint=NA, legend.show=F) +
  tm_layout(frame = F)
tmap_save(map, "map.png")

## other visual works over adobe


###############################################################################
############################ Sensitivity Analysis #############################
###############################################################################

pacman::p_load(sensemakr)

summary(lm2)
lm2_sens_1 <- sensemakr(model=lm2,
                        treatment="treatB",
                        benchmark_covariates=c("q2_polbias","q3_gap"),
                        kd=1:3)

lm2_sens_2 <- sensemakr(model=lm2,
                        treatment="treatC",
                        benchmark_covariates=c("q2_polbias","q3_gap"),
                        kd=1:3)

summary(lm2_sens_1)
summary(lm2_sens_2)
